﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Mclib
{
	/// <summary>
	/// A class that approximates calling Math.Exp() but evaluates the function more quickly.  Precision is imperfect
	/// as this requires using linear interpolation and a lookup table.
	/// </summary>
	public class FastExp
	{
		public double[] LookupPos_0_1 = new double[1001];
		public double[] LookupPos_1_10 = new double[901];
		public double[] LookupPos_10_100 = new double[901];
		public double[] LookupPos_100_709 = new double[610];
		public double[] LookupNeg_0_1 = new double[1001];
		public double[] LookupNeg_1_10 = new double[901];
		public double[] LookupNeg_10_100 = new double[901];
		public double[] LookupNeg_100_709 = new double[610];
		/// <summary>
		/// Constructs a new instance of the FastExp class.  This should be called as few times as possible because it requires
		/// approx 3500 evaluations of the Math.Exp() function and approx 3500 division operations.
		/// </summary>
		public FastExp()
		{
			int i;
			double z, ez;
			for(i=0; i<1001; i++)
			{
				z = (double)i * 0.001;
				ez = Math.Exp(z);
				LookupPos_0_1[i] = ez;
				LookupNeg_0_1[i] = 1.0/ez;
			}
			for(i=0; i<901; i++)
			{
				z = (double)(i+100) * 0.01;
				ez = Math.Exp(z);
				LookupPos_1_10[i] = ez;
				LookupNeg_1_10[i] = 1.0/ez;

				z = (double)(i+100) * 0.1;
				ez = Math.Exp(z);
				LookupPos_10_100[i] = ez;
				LookupNeg_10_100[i] = 1.0/ez;
			}
			for(i=0; i<610; i++)
			{
				z = (double)(i+100);
				ez = Math.Exp(z);
				LookupPos_100_709[i] = ez;
				LookupNeg_100_709[i] = 1.0/ez;
			}
		}
		/// <summary>
		/// Evaluates Exp(z) by using linear interpolation through lookup tables.  This is faster than directly evaluating
		/// Math.Exp(z) but not as precise.
		/// </summary>
		/// <param name="z">The operand in (-Inf, +Inf).</param>
		/// <returns>Exp(z).</returns>
		public double Eval(double z)
		{
			double wHigh;
			int iLow;
			if(z<0.0)
			{
				if(z>-10.0)
				{
					if(z>-1.0)
					{
						// z is [0, 1.0)
						z*=-1000.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupNeg_0_1[iLow]*(1.0-wHigh) + LookupNeg_0_1[iLow+1]*wHigh;
					}
					else
					{
						// z is [1.0, 10.0)
						z = (-z-1.0)*100.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupNeg_1_10[iLow]*(1.0-wHigh) + LookupNeg_1_10[iLow+1]*wHigh;
					}
				}
				else
				{
					if(z<=-709.0)
						return 0.0;
					if(z>-100.0)
					{
						// z is [10.0, 100.0)
						z = (-z-10.0)*10.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupNeg_10_100[iLow]*(1.0-wHigh) + LookupNeg_10_100[iLow+1]*wHigh;
					}
					else
					{
						//	z is [100.0, 709.0)
						z = -z-100.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupNeg_100_709[iLow]*(1.0-wHigh) + LookupNeg_100_709[iLow+1]*wHigh;
					}
				}
			}
			else if(z>0.0)
			{
				if(z<10.0)
				{
					if(z<1.0)
					{
						// z is [0, 1.0)
						z*=1000.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupPos_0_1[iLow]*(1.0-wHigh) + LookupPos_0_1[iLow+1]*wHigh;
					}
					else
					{
						// z is [1.0, 10.0)
						z = (z-1.0)*100.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupPos_1_10[iLow]*(1.0-wHigh) + LookupPos_1_10[iLow+1]*wHigh;
					}
				}
				else
				{
					if(z>=709.0)
						return double.PositiveInfinity;
					if(z<100.0)
					{
						// z is [10.0, 100.0)
						z = (z-10.0)*10.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupPos_10_100[iLow]*(1.0-wHigh) + LookupPos_10_100[iLow+1]*wHigh;
					}
					else
					{
						//	z is [100.0, 709.0)
						z -= 100.0;
						iLow = (int)z;
						wHigh = z % 1.0;
						return LookupPos_100_709[iLow]*(1.0-wHigh) + LookupPos_100_709[iLow+1]*wHigh;
					}
				}
			}
			return 1.0;
		}
	}
}